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Phase Field Crystal Study of Symmetric Tilt Grain Boundaries of Iron 

A. Jaatineni'2, C. V. Achim^, K. R. Elder*, T. Ala-Nissilai^^ 



We apply the phase field crystal model to study the structure and energy of symmetric tilt grain boundaries ofbcc 
iron in 3D. The parameters for the model are obtained by using a recently developed eight-order fitting scheme [A. 
Jaatinen et al., Phys. Rev. B 80, 031602 (2009)]. The grain boundary free energies we obtain from the model are 
in good agreement with previous results from molecular dynamics simulations and experiments. 



1 Introduction 



One of the challenges in modern computational materials science is being able to access phenomena that take place 
on different time and length scales. Recently, a new model called the phase field crystal (PFC) mo del has been 
const ructed that describes phenomena taking place on atomic length scale but diffusive time scales jElder et al 



20021) . the combination of which has been unaccessible for molecular dynamics (MD) simulations using present 
day computers. The PFC method is able to achieve this combination of scales by replacing the individual atoms in 
MD simulations by a continuous field that exhibits the periodic nature of the underlying atomic lattice in the solid 
phase, but evolves in diffusive time scales. As the original formulation of the PFC model relies on a phenomeno- 
logically constructed free energy functional, it has proven challenging to quantitatively relate the parameters that 
enter the PFC model to the properties of real materials. In a recent study, we have shown how the PFC model can 
be modified in such a way that many relevant material properties can be included in the PFC free energy, while still 
conserving most of the computational simplicity of the model (Jaatinen et al., 2009). This can be achieved using 
the so-called Eighth Order fitting (EOF) scheme. In the present work, we use the EOF-PFC model to study the 
energy of symmetrically tilted grain boundaries of bcc iron near its melting point. Results of our calculations will 
provide an independent theoretical prediction for the grain boundary energy of iron near its melting point. Perhaps 
more importantly, they also provide indication on whether the EOF sche me is a good candidate for future studies of 
grain boundary related phenomena, such as grain boundary premelting ("Berrv etaL , 2008 ; Mellenthin et al. , 2008 ) 
and nanocrystalHne deformation (.Elder and Grant..2004:,Stefanovic et al...2009l) . with the PFC model. 



2 Model 



The phase field crystal (PFC) model describes crystallization phenomena and has been applied to many differ- 
ent problems in materials s cience. Originally, the model was phenomenologically postulated in the spirit of the 
Ginzburg-Landau theory bv lElder et al. f (l2002h . The model consists of an order parameter field that is driven by 
dissipative, conserved dynamics to minimize a free energy functional, whose ground state exhibits a periodic struc- 
ture commensurate with the crystal symmetry of interest. More specifically, the deterministic equation of motion 
for the locally conserved order parameter field n(r) in these models is given by 
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where M is mobility and F is free energy of the system as a function al of the field n(f, t). The mo st usual choice 



of a free energy functional for PFC studies is the form derived by ISwift and Hohenberg (119771) for a study of 

for obtaining a 
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convective instabilities. Here we briefly describe the procedure suggested by 
free energy that is capable of describing body-centered cubic materials from a fundamental basis. 



D erivation of our free energy starts fro m the classical density-functional theory (DFT) of freezing, first pioneered 
bv lRamakrishnan and Yussoufn ( 1197 9). The simplest free energy used in DFT of freezing studies can be obtained 
by expanding the excess (over ideal gas) contribution to the free energy around a uniform reference density po. 



leading to 
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where Ai^ is the free energy difference between a given state and the reference, ksT is the thermal energy and 
the order-parameter field n is related to the ensemble-averaged one-particle density p{r) through n{r) — (p(r) — 
Po)/Pq. The function C entering Eq. (|2]i is the two-body direct correlation function (DCF) of the reference liquid, 
which is related to the total pair correlation function, and thus structure factor of the liquid through the well-known 
Ornstein-Zernike relation. 

In principle it is possible to use Eq. (|2]l for PFC studies; however, in practice solutions of this equation are very 
sharply peaked in space (around atomic lattice positions) which restricts numerical calculations to very small 
systems. To overcome this difficulty, we make two approximations. Firstly, the non-local part of F is simplified 
by expanding the DCF in fc-space up to eighth order as 

dk) « c{k,.) r (^^)' - Es (^%^)', (3) 

where the parameter fc„j is the position of the first maximum in the original DCF, and F and Eb are chosen such 
that the expansion reproduces the position, height and curvature of the first maximum, and the fc = value of the 
original DCF. This is achieved by choosing 
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where primes denote derivatives with respect to k, and 

Eb = C{Kn) ~ C(0) - P. (5) 



As shown by Ijaatinen et al. ( 20091) . this eighth order expansion provides an excellent fit from fc = to the first 
peak in C{k) in the case of iron. However, for the density oscillations whose wave vector corresponds to a fc larger 
than km, the expanded free energy now results in a significantly larger free energy penalty, thus reducing stability 
of the solid phase. In order to partly correct this error, and to further simplify the mathematical form of the free 
energy, the local part is expanded in a fourth order power series as 
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The coefficient 1/2 for the second order term in this series is obtained from a Taylor series in order to conserve 
the linear stability of the model. For the third and fourth order terms we have included constants a and b, which 
are chosen such that the solid will stabilize in the model, with the correct amplitude of density fluctuations cor- 
responding to the first star of reciprocal lattice vectors, Ug. This is approximately achieved by choosing them 
as 

(7) 
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where iS'(fc„i) = (1— C(fcm)) ^ is height of the first peak in the structure factor of the reference liquid dJaatinen et al 
E009.) . With these approximations, the free energy finally becomes 
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In our previous study we fitted the parameters of the current mod el with molecular dynamic s si mulation data from 
Wu and Karma ( 2007 ). complemented by experimental data from Itami and Shimoii ( 1984 ) and Jimbo and CrambI 
(Il993b . for iron. The resulting model parameters are fc,„ = 2.985 A-\ C{km) = 0.668, T = 11.583 and Eb = 
38.085 for the DCF part, and a = 0.6917 and b — 0.08540 for the non-linear part. The reference temperature and 
density, necessary for converting model results into SI units are T = 1772 K and po — 0.0801 A^'^. We have also 
shown that with these parameters, the model produce s very reasonable pre dictions for bulk modulus, solid-liquid 
coexistence gap and low index surface free energies. (!Jaatinen et al. . 20091) In the next section, we will show that 
this model is also capable of describing grain boundaries in iron. 




Figure 1: Grain boundary free energy as a function of the misorientation angle 9 when the tilt axis is in the (100) 
direction. Solid line shows the best fit to the Read-Shockley equation at small tilt angles. Dashed line is a guide to 
the eye. 



3 Grain boundaries 



Using the model and parameters described in the previous section, we have studied the symmetrically tilted grain 
boundaries of iron. Three choices of tilting axis, (100), (110) and (111), were studied. The calculations were 
performed at a density of [p — pQ)/pa = 0.030, which is just above the density of the solid at coexistence. 
A density close to the solid-liquid coexistence was chosen due to the fact th at the model parame ters have been 



obtained to reproduce the selected properties of iron close to its melting point (iJaatinen et al.l 12009 ) 



First, we have found the lattice constant a of the bcc phase that minimizes the free energy at the density in question. 
Then, a computational box of dimensions [L^, Ly, L^), and periodic boundary conditions in all directions was 
initialized by using the one-mode approximation, 

n{r) = no + {1 + no)4u [cos{qx) cos(gy) + cos(qx) cos(gz) + cos(gj/) cos{qz)] , (9) 

where q = 27r/abcc, Obcc being the lattice spacing, and u is the amplitude of density fluctuations. First, a rotation 
to the one-mode approximation was applied such that the tilt axis was parallel to the z axis of the computational 
box. Then, in the region a; = . . . Lx/2, the one-mode approximation was rotated in the {x, y) plane by an angle 
—9/2, and in the region x = Lx/2 . . . hy an angle 9/2, creating two symmetrically tilted grain boundaries, 
located at .t = and x — Lx/2, with a tilting angle 9. The tilting angles 9 and box dimensions {Lx,Ly, Lz) were 
chosen such that the periodic pattern of the lattice continues smoothly over the periodic boundaries in the y and 
z directions, and an integer number of atomic planes was placed between the two grain boundaries in the system. 
The system size in the direction of the tilting axis, L^, was chosen to be the closest distance between two similar 
atomic planes in the z-direction. In the perpendicular directions, and Ly were varied from case to case, such 
that for small angles, the simulation box contained several thousands of maxima, while for large angle calculations, 
only hundreds of maxima were required, due to the smaller spacing between dislocations at the boundary. 

In order to find the minimum energy of the system with the grain boundaries , we have numeric ally integrated 



Eq. ([T]l until convergence, using the semi-implicit operator splitting method of iTegze et al.l (l2009h . The Laplace 
operator is discretized in fc-space as Ak = —k^. Grain boundary free energy per unit area, Egb, is then obtained 
from 

Egb = Lz^-^, (10) 

where / is the average free energy density of the final configuration, and is free energy density of a non-rotated 
system. 



Fig. [T]shows the results of our grain boundary free energy calculations, and Fig. |2]shows three examp le final con- 
figurations for the (100) case. As in previous studies utilizing the original PFC free energy functional dElder et al. , 
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Figure 2: Example final configurations from the (100) case. The gray scale in these figures corresponds to 
n{x,y,0). On the left, we show the low-angle bcc(0 64 1) (100) i;40975 (9 = 1.79°) boundai-y, where one 
clearly sees the individual dislocation. The middle figure shows a large angle bcc(051) (100) S13 (6 = 22.62°) 
boundary, and the figure on the right shows the bcc(031) (100) E5 (6 = 36.87°), for which a small energy cusp 
was observed. 
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Figure 3: Grain boundary free energy as a function of the misorientation angle 6 when the tilt axis is in the (110) 
direction. Dashed line is a guide to the eye. 
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Figure 4: Example final configurations from the (110) case, with the gray scale corresponding to n{x, y, 0). From 
left to right, the corresponding grain boundaries are bcc(32 32 1) (110) i;2G49 (9 = 2.53°), bcc(331) (110) E19 
(9 = 26.52°), bcc(332) (110) Sll (9 = 50.47°) and bcc(112) (110) E3 (9 = 109.47°). 



2002t Elder and Grand l2004i iMellenthin et al.l l2008h. we find that for small angles the grain boundary energy 
closely follows the expression by Read and Shocklevn 1950l) . 



Egt = Eoe{A~\n (0)), 



(11) 



where Eq is an energy scale that depends on elastic properties of the material studied, and the constant A depends 
on the dislocation core energy. In addition, both of these constants depend on the orientation of the grain boundary. 
In the present work, as w e have not developed num erical tools for quantifying all the properties that enter the 
expressions for Eq and A ( Read and Shocklev , ll950l) . we have chosen to treat them as fitting parameters. The best 
fit to the current data gives Eq = 1.1271 Jm"^ and A = -0.41. 

For large values of 9, the changes in Egi, are moderate, with the maximum value being 0.377 Jm~^. By close 
inspection of Fig. [T]one can distinguish small energy cusps at angles corresponding to bcc(031) (100) S5 (6 = 
36.87°) and bcc(053) (100) S17 (9 — 61.93°) grain boundary planes. Both of these energy cusps correspond 
to a low coincidence site lattice (CSL) E. However, for the other low-S grain boundaries, for example the 
bcc(021) (100) E5 (9 = 53.13°), no such cusps were observed. More generally, we do not observe any clear 
corr elation between the CSL S and energies of large angle grain boundaries. Thus, our results do not agree 
with IZhang et alj (2005) who found, using molecular statics with modified analytical embedded atom method, 
that the energy of symmetrically tilted grain boundaries in iron increases with increasing S. Our results agree 
much better with Shibuta et al. (2008) who have performed molecular dynamics simulations with Finnis-Sinclair 
potential to calculate the energies of different symmetric tilt boundaries of iron in different temperatures. Their 
results are qualitatively (but not quantitatively) very similar to ours, and they also observe the small ene rgy cusp 
at bcc(031) (100) S5 present in our results as well. However, unlike in our results, [shibuta et aL I (l2008h observe 
another small energy cusp at bcc(021) (100) S5. The cusp we observe at bcc(053) (100) S17 is not present in their 
results. 



Results of our grain boundary free energy calculations for the (110) case are shown in Fig. [5] and example 
final configurations for this case are shown in Fig. |4] The low-angle behavior of Egb vs. 9 is qualitatively 
very similar to the (100) case, but at large angles, interesting features are observed. There is a deep minimum 
corresponding to bcc(112) (110) S3 {9 ~ 109.47°), which is a twin boundary. The grain boundary energy in this 
minimum is only approximately third of the maximum values observed. Other clear local minima are observed at 
bcc(221) (110) S9 (9 = 38.94°), bcc(lll) (110) S3 (9 = 70.53°) and bcc(114) (110) S9 (9 = 141.06°). Thus 
for (110) grain boundaries all angles with S < 10 correspond to a local minimum in grain boundary energy. 
However, as in the (100) case, we do not find a more general correlation between S and j^oj,, as can b e seen from 
Fig|5] where we plot Egi, vs. S for the (110) case. In comparison of our results with IShibuta et alj ( l2008h . we 
find that our depth of the bcc(112) (110) S3 minimum is in an excellent agreement with their work. However, 
their bcc(221) (110) S9 minimum is much less pronounced than in our work, and the remaining two of the local 
minima, that we observed, were not observed in their study at all. Instead, they observe a small energy cusp at 
bcc(554) (110) S33 (9 = 58.99°), which we did not find anomalous. 

Results for the case where tilting axis is in the (111) direction are shown in Fig. |6] In this case, the energy 
increases with increasing misorientation, without any noteworthy energy cusps, until reaching a maximum value 
at 9 = 32.20°. When further increasing 9 from that maximum, the grain boundary decreases, as the orientation 
gets closer to the partial matching at 6* = 60°. Finally, at 6* — 60° (i.e. bcc(121) (111) E3) the grain boundary 
energy is 0.357 Jm~^, which is approximately 30 % of the maximum value. The place of the maximum energy, 
magnitude of the 9 — 60° value c ompared to the maxi mum, and the absence of energy cusps agree very well with 
the molecular dynamics results of I Shibuta et al.l (120081) . 



Exper imentally, the large angle grain boundary energy of S iron has been determined to be 0.468 Jm^^ dMun , 
19751). in good ag reement with our maximum value. It is also worthwhile noting that the agreement of our results 
with lMurn (119751) for the ratio of large angle grain boundary energy to the solid-liq uid surface free ene rgy is even 
better: Enh/cTs -i ~ 2.16 in the present model (solid-liquid surface energy from Ijaatinen et aP (Eob9) ). where 
the values from iMun ( 1975 ) give Egb/og-i ~ 2.29. When comparing our values for grain boundary energies 



with the results of the pr eviously mentioned atomistic calculations, we observe that the grain boundary energies 
obtained bv Zhang et al. 1 (12005) are of the order of several Joule s per sq uare meter, being significantly higher than 



th e values reported h ere or those reported by by IShibuta et al.1 (l2008h . It should be noted that the calculations 



of IZhang et al. I (12005 1) were conducted at absolute zero temperature, no t near th e melting point, which perhaps 
explains this discrepancy. The grain boundary energies reported by Shibuta et al.l (|2008 ). on the other hand, reach 
a maximum of approximately 1.6 Jm^^ at moderate temperatures, rising to approximately 2.4 Jm^^ close to the 
melting point. A likely reason for our grain boundary energies being closer to the experimental values than those 
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Figure 5: Grain boundary free energy as a function of the coincidence site lattice S when the tilt axis is in the 
(110) direction. 
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Figure 6: Grain boundary free energy as a function of the misorientation angle 9 when the tilt axis is in the (111) 
direction. Dashed line is a guide to the eye. 



of Shibuta et al. ( 2008h is that we are, in accordance to the experimentaHsts, actually measuring the grain boundary 
free energy, unlike 'S hibuta et aL (2008) who only measure the excess energy of creating a grain boundary, without 
taking into account the entropic part of the free energy. 



4 Conclusion 

We have applied the recently established EOF-PFC model to study the structure and energy of symmetric tilt 
grain boundaries in iron, with the axis of rotation parallel to (100) and (110) directions. For small misorientation 
angles, we found that the grain boundary energy increases with increasing misorientation, in agreement with the 
Read-Shockley equation. For large misorientation angles, we observed local minima in the grain boundary energy 
vs. misorientation curve, many of which have been observed in molecular dynamics simulations as well. The 
local minima are associated with low values of the coincidence site lattice S, but no general correlation between 
Yj and the grain boundary energy was found. In general, the qualitative agreement between present results and 
previous MD results at finite temperature was found to be good. Quantitatively, the agreement of our result with an 
experimental value for the large angle grain boundary energy is better than that from the MD. As our phase field 
crystal results agree well with previous MD study of grain boundaries in iron, and an experimental large angle 
grain boundary energy, they may even be considered an independent prediction for the angular dependence of 
grain boundary energy in iron. More importantly, they show the potential of the EOF-PFC in modeling materials 
properties where grain boundaries play an important role. 
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